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CONSTRAINED BURN OPTIMIZATION FOR THE 
INTERNATIONAL SPACE STATION 


Aaron J. Brown* and Brandon A. Jones’ 


In long-term trajectory planning for the International Space Station (ISS), trans- 
lational burns are currently targeted sequentially to meet the immediate trajec- 
tory constraints, rather than simultaneously to meet all constraints, do not employ 
gradient-based search techniques, and are not optimized for a minimum total delta- 
v (Av) solution. An analytic formulation of the constraint gradients is developed 
and used in an optimization solver to overcome these obstacles. Two trajectory ex- 
amples are explored, highlighting the advantage of the proposed method over the 
current approach, as well as the potential Av and propellant savings in the event 
of propellant shortages. 


INTRODUCTION 


The first module of the ISS, Zarya (“sunrise” in Russian), was launched on November 20, 1998 
from Baikonur Cosmodrome in Kazakhstan. It was launched into a 400 km altitude, roughly circular 
orbit, and an inclination of 51.6°. Each module that followed was launched into the same orbit and 
added to the ISS main body. With only minor fluctuations in altitude (between roughly 330 km and 
420 km), the ISS has ostensibly remained in the same orbit ever since. 


The primary perturbing forces acting on the ISS orbit include Earth’s non-spherical gravity, third- 
body perturbations from the Sun and Moon, and atmospheric drag. Except for the secular effect of 
Jz and atmospheric drag, the effects of these forces on the ISS orbit are either relatively small in 
magnitude, periodic in nature, or both, and therefore do not affect the orbit’s long-term evolution. 


The ISS orbit is planned and maintained in part by the ISS Trajectory Operations (TOPO) group 
at the NASA Johnson Space Center. TOPO personnel routinely generate look-ahead trajectory plans 
for the ISS, ranging in length from six weeks to almost two years. Each plan contains the trajectory- 
related events that occur during the plan’s time frame, including visiting vehicle launches, dockings, 
undockings, landings, and ISS translation burns. 


One of the goals in generating a given plan is to target the ISS translation burns that satisfy the ISS 
trajectory constraints. The burns are scheduled to occur on agreed-to dates by the ISS International 
Partners (IPs). The constraints are imposed by Russian Soyuz (crewed) and Progress (cargo) vehicle 
launches and Soyuz vehicle landings, which require the ISS orbit to fall within a specified Longitude 
of Ascending Node (LAN) range on the day of the launch or landing. The ISS LAN itself shifts 
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approximately 23.3° westward per orbit due to the combination of Earth rotating to the east and 
the ISS orbit precessing to the west. Additionally, one or more semi-major axis (SMA) altitude 
constraints may be imposed on the trajectory. These constraints ensure the ISS achieves a certain 
mean SMA altitude on a certain date, or maintains a mean SMA altitude across a range of dates. 


Unfortunately, the legacy software used by TOPO personnel for targeting these burns does not 
utilize any gradient-based search techniques for finding feasible burn solutions. This makes for a 
cumbersome process in which burns are targeted manually rather than automatically, thereby forcing 
users to find burn solutions sequentially (1.e, burn by burn), rather than simultaneously for all burns. 
Furthermore, this process does not at all enable finding optimal solutions that minimize total Av 
over the trajectory plan time frame. 


In recent years, however, a large-scale project has been underway to replace this legacy software 
with a set of trajectory tools that can utilize differential correction to find feasible burn solutions, 
or can be embedded in a larger optimization framework to find optimal burn solutions. The goal 
of this work is to provide the mathematical formulation of the constraint gradients in the ISS burn 
optimization problem, and to use these gradients in an optimization solver to find either feasible or 
minimum Av burn solutions for ISS trajectory planning. 


Reboost Burn Modeling 


In ISS trajectory planning, translation burns are modeled in the Local- Vertical, Local Horizontal 
(LVLH) reference frame. The LVLH frame is centered on the vehicle, with the Z-axis (the LV axis) 
opposite the radius vector (i.e. down), the Y-axis opposite the orbit’s angular momentum vector, 
and the X-axis (1.e. the LH-axis) completing the right-handed system, as shown in Figure 1. 
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Figure 1: LVLH Coordinate System! 


Each burn is modeled as impulsive, with Av > 0, and in a fixed direction defined by LVLH yaw 
(w) and pitch (@) angles in a yaw-pitch-roll sequence. These angles are typically small, resulting in 
a posigrade or reboost burn. Reboost burns serve to restore orbit altitude lost due to atmospheric 
drag. Each burn is scheduled to occur on a specific date, with the burn’s Time of Ignition (TIG) on 
that date occurring at apogee (v = 180°) of Daily Orbit 1 (DO1). Since the ISS orbital period is 
roughly 92 minutes, there are either 15 or 16 DOs each day. On a given day, DO1 is defined to be 
the first orbit with a Longitude of Ascending Node (LAN) that is west of 20°E. Burns are nominally 
scheduled to occur on DOI to allow for their visibility via Russian ground stations. 


Burn Optimization Problem Statement 


Given the ISS reboost burn model in the previous section, the burn optimization problem can be 
stated as follows. Given a fixed timespan from to to t, with mn reboost burns to occur on fixed dates 
Gay eee 

mr 
minimize J = >». Avj. (1) 
i=1 
Since each Av; > 0, J in Equation | is effectively the @;-norm of an n-element vector containing 


the Av magnitudes. 


The controls in this problem are the burn Av magnitudes and the constraints, as noted in the 
Introduction, are given by the Russian Soyuz and Progress launches, and Soyuz vehicle landings. 


e Soyuz launch LAN constraints: 13.9° < A < 15.4° 

e Soyuz landing LAN constraints: 12.7° < A < 16.5° 
e Progress launch LAN constraints: 26.5° < A < 38.9° 
e SMA altitude constraints 


Given the cost function J in Equation 1, the objective partial derivatives are easily given by 


OJ 
——_—=], i=l,...,n. Z 
The constraint partial derivatives are now developed by first examining the total state differential for 
a generic trajectory segment. This differential is used repeatedly to formulate state partial deriva- 


tives, which are then used to arrive at the desired constraint partial derivative expressions. 


Total State Differential for a Generic Trajectory Segment 
The development in this sub-section follows and expands on work done by Ocampo and Munoz.? 


Consider the generic trajectory segment in Figure 2. Let the state x be defined as 
x2[rF vw ]’, (3) 


where r and v are the J2000 position and velocity vectors. Furthermore let dx; be the total state 
differential at t¢ following Ax,y. The goal of this sub-section is to express dx; as a function of 


the independent variables dtp, dx , d(Axo), dtr, and d(Ax,). The final result for dx; is used 
repeatedly in the next section to develop state partial derivatives for a notional set of ISS trajectory 
events. 


Begin by noting that the time-fixed variations 6 Xp and 3 xp are related through the state transition 
matrix (STM) ®(t;+, to), abbreviated as ®¢ 9, through 


OXp = Dro 6x¢. (4) 


Next note that the plus and minus states at fp and ty are related through Axo and Ax,, and also 
lead to the following total differential expressions for dx¢ and dx}. 


a =x) +Axp => dx¢ = dx, + d(Axo) (5) 
x =x, + Ax —s dx; = dx, +d(Axy) (6) 
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Figure 2: Generic Trajectory Segment, as depicted by Ocampo and Munoz.” The state x evolves 
as: (to, xp) + Axo + (to, xd) > (tp, x7) > Axy > (ty,xf ] 


The total differentials dx, and dx¢ are also related to their respective time-fixed variations through 


dx, = OX + x dts (7) 


dx) =6xj +xfpdtop => Oxf =dxf — xf dto. (8) 
Using Equations 4 through 8, it is readily shown that 
dx; = fo |dxq — xydto + d(Axo)| + x;dt; + d(Axy). (9) 


which is the desired result in terms of the independent variables. Ultimately though, the real inde- 
pendent variables, or controls, in this problem are the burn Av magnitudes as noted previously. The 
quantities Axo and Ax; in Equation 9 therefore become dependent variables, themselves functions 
of the burn Av magnitudes. This effect is addressed in the next section on state partial derivatives. 


State Partial Derivatives 


The state partial derivatives are best related by examining a notional set of ISS trajectory events, 
as shown in Figure 3. This figure depicts seven trajectory events that occur at times ¢; through t7, 
with an initial coast arc that occurs from tg to t;. Reboost burns (Av) occur at Events 1, 3, and 6, 
while LAN constraints (A) occur at Events 2, 4, 5, and 7. 


State differentials are now developed, event by event, using Equation 9 as a guide. The state 
partial derivatives are obtained directly from these expressions and used in the next section to form 
the desired constraint partial derivatives. 


Event 1 


Applying Equation 9 to Event 1, and assuming xo, to, and ¢; are fixed, thendx) = d(Axo) = 0, 
and dtp = dt; = O, leaving 
dx} = d(Ax}). (10) 


However, as noted in earlier, Ax, is itself a function of Av;, namely 
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Ke a | (11) 
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Figure 3: Altitude (/) vs. time (t) for a notional set of seven ISS trajectory events. These events 
include reboost burns (Av) and LAN constraints (A) that must be met by those burns. 


where fue is the transformation matrix from the LVLH frame (L) to the J2000 frame (J), and ti is 
the Av unit vector in the LVLH frame defined by 7 and 91, 


a, = cos 6; cosw, cos@;sinyw, sin@, ie (12) 


Equation 10 therefore becomes 


ix; = _ | d(Av1). (13) 


Event 2 


Applying Equation 9 now to Event 2, with x, and ¢, fixed, and Ax. = 0, 
dx 9 = O21 d (Ax) + X9 dto. (14) 


Ignoring the x2 dtz term for a moment, the rest of Equation 14 takes an arbitrary perturbation in 
Ax}, i.e. d(Ax;), and maps it via the STM to the corresponding perturbation in xg, i.e., dxe. 
Now along the nominal, unperturbed trajectory, x2 occurs at an ascending node, since the LAN 
constraint at Event 2, by definition, requires evaluating the spacecraft’s longitude when it is at the 
ascending node. Any perturbed trajectory must likewise have its x2 occurring at the ascending 
node in order to facilitate an apples-to-apples comparison of LAN values between the perturbed and 
nominal trajectories. 


Now being at the ascending node in both the nominal and perturbed cases requires the z-component 
of ra (1.e. Z2) to be 0 in the True-of-Date (TOD) Earth Equator coordinate system. The TOD system 
is similar to J2000, but is defined using the true Earth equator and true equinox of date. While the 
ascending node can be referenced to either J2000 or TOD, the TOD system is used since the LAN 
by definition is measured along Earth’s true equator. The z2 = O requirement is already satisfied 
by definition along the nominal trajectory. Along the perturbed trajectory, however, d (Ax ,) prop- 
agated to the same fz results in a perturbed x2 that is no longer at the ascending node, falling above 
or below the node, depending on the nature of d(Ax,). Left uncorrected, the perturbed x2 cannot 


be compared with the nominal trajectory since it is not at the ascending node, thus leaving the LAN 
undefined, or at least unknown. 


To remedy the situation, the x2 dt2 term is re-introduced and utilized to enforce the z2 = O 
requirement, thereby bringing the perturbed x2 up to or down to the ascending node as needed. 
This 1s accomplished by first expressing the z2 = 0 requirement more formally as 


T 0,, 
z=kix.=0, kL 4[0 0100 0) oe ee (15) 
03,3 ‘a 


where de is the 3x3 transformation matrix from J2000 to TOD, the A subscript denotes the 
z2 = O requirement on LAN, and the event number (2) sub-subscript defines the epoch at which 


TOD , ‘ 
T jor 1S evaluated. Proceeding, 


29 =ky,x2=0 => dz =kj,dx2 =0, (16) 


and Equation 14 therefore becomes 


k\,dx2 = 0 = kj, ®21 d(Ax)) + kj, x2 dtp. (17) 
Solving Equation 17 for dt, 


T 
k\, 
ae e 
ky x2 


dtz = — 


O21 d(Ax}). (18) 


dx» 1s then found by substituting dt2 back into Equation 14 to obtain 


: T 
x2k) , 


O21 d(Ax}). (19) 


dx2 = i ee ; 
1 : kj, x2 


The preceding development is best understood by re-examining Equation 14 in light of the time- 
fixed variation 0x, 


dxg = 0X2 + xX9dtg, 6x2 = P21 d (Ax;), (20) 


and in conjunction with Figure 4. This figure shows a nominal trajectory that stops at the ascending 
node and a perturbed trajectory that stops above the ascending node, both at time tg. The quantity 
0X2 represents the time-fixed variation in x2. The quantity dtz then represents a time slip in fg, 
computed such that dxyg = 0x2 + x2 dtg lies along the TOD equator, thus satisfying the dzz = 0 
requirement. This process is dubbed “ascending node shaping,” since it takes an x9 that 1s initially 
dispersed in all six components, and flattens the dispersion such that dx lies along the TOD equator. 
As mentioned earlier, this shaping is necessary in order to properly compare LAN values between 
the perturbed and nominal trajectories. 


Defining the following shaping matrices, 


ki Xok} 
he = Ee ; and Sh, < 1 7 Ee — ) (21) 
do 2 do 2 
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Figure 4: Ascending Node Shaping 


Equations 18 and 19 become 


dtz = S} & 1 d (Ax) (22) 
dx2 = S22) d (Ax). (23) 
The terms “time slip” and “shaping matrix” are borrowed from Moesser,’ who uses Equations 18 


and 19 in the context of linear covariance analysis to examine the Lunar powered descent problem. 
Moesser further notes that Ss. is an idempotent shaping matrix, which exhibits the property 


(s3,) =5),5° =). (24) 


Another property worth noting is that Se. contains at least one zero eigenvalue and is therefore 
rank deficient. This can be deduced geometrically by noting that S 0 , always zeros out (and therefore 
removes information from) the z-channel by enforcing z2 = 0. 


Finally, substituting Ax, from Equation 11 into Equations 22 and 23 yields 


0 
dty = S) bo _ d (Av) (25) 
T7(r1,V,) Uy 
d O3x1 
ax9 = Sx, P21 d (Av}). (26) 


Ti (r1,v,)- ty 
These expressions can also be written more compactly by utilizing 6x in Equation 9 again to yield 
dtz = Si. dx2 (27) 
dx = Ss. 0X9. (28) 


This form of the time and state differential expressions is useful since it is identical for subsequent 
events and therefore reduces code complexity. 


Event 3 


Applying Equation 9 now to Event 3, with Axy = 0, 
dx} = ®3 9 [dxg — X2 dtg] + x3 dt3 +d (Axs). (29) 
Substituting dt2 from Equation 18 and dx2 from Equation 19 into Equation 29 yields 


dxf = 39 ®21d(Ax,) + x3 dt3 + d(Axs). (30) 


At this point a similar shaping analysis must be performed for Event 3. This event, a reboost 
burn, does not occur at an ascending node as in Event 2, but rather at apogee, and is therefore 
dubbed “apogee shaping.” Whereas the requirement by ascending node shaping was z2 = O in the 
TOD system, the analogous requirement by apogee shaping is 


COs ’y3. = 13 Vz = 0. (31) 
This requirement, and the shaping derived from it, ensures that the reboost burn is in fact performed 
at apogee. Left uncorrected, a perturbed x, resulting from d(Av;) would not be at apogee, and 


therefore the reboost burn would not occur at apogee as needed. Mathematically, this requirement 
also holds at perigee. This ambiguity, however, is eschewed by noting that 


O (rove _ = 
ONS) o 3) =r3Vv, +(v3)° (32) 


is negative at apogee and positive at perigee. Proceeding, 
(casa, ) =k dx, 0, ke) ae |, (33) 
where the +y subscript denotes the cos y requirement at apogee. Returning to Equation 30, 
dx} = dx; +d(Axs), dxz = 32 21 d(Ax;) + x3 dts. (34) 
Therefore, 


k) dx; =0 =k}, 032 21 d(Ax)) +k), x; dts. (35) 


3 3 


Following a procedure similar to Event 2 results in 


kf 
dt3 = aT ae O39 O21 d (Ax; ) (36) 
73 7*3 
x3 k;, 
dX = | Pee = KT a O39 O21 d(Ax}). (37) 
3 7"3 


As before, the preceding development is best understood by examining dx, in Equation 34 in light 
of 0x, 
dX = 0X3 ++ Xe dts, 0X3 = O39 O21 d (Ax;), (38) 


and Figures 5 through 7. Figure 5 shows a nominal trajectory that stops at apogee and a perturbed 
trajectory that considers only position dispersions about apogee, both at time t3*. The quantity ors 
represents the time-fixed variation in rz, and r3 represents r3 along the perturbed trajectory. The 
quantity dt3 then represents a time slip in t3, computed such that drg3 = dr3 + rz dt3 is orthogonal 


to v, , thus ensuring that (v5 ) "drs = 0 is Satisfied. Figure 6 shows the opposite case, considering 
only pre-burn velocity dispersions at apogee. dts in this case is the time slip in ¢3 required such that 
dV; = 0V3 + Vz dt3 is orthogonal to rz, thus ensuring that (r3) dv = (is satisfied. Finally, 
Figure 7 shows both position and pre-burn velocity dispersions at t3 resulting from d(Av,). dts 
is now the time slip in ¢3 required to ensure d (r$ v5 ) = 0 1s satisfied. This enables the burn at 
Event 3 to be performed at apogee along the perturbed trajectory (denoted by (r3),, and (v5 ). in 
Figure 7), as desired. 
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Figure 5: Apogee Shaping, Considering Position Dispersions 
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Figure 6: Apogee Shaping, Considering Pre-Burn Velocity Dispersions 
“Note that this is a fictitious perturbed trajectory, since a perturbation in Av, d (Av), would result in both position 


and pre-burn velocity dispersions at t3, as shown in Figure 7. Considering these two dispersions separately, however, 
aids in understanding the mechanism of apogee shaping in the combined dispersion case. 
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Figure 7: Apogee Shaping, Considering Position and Pre-Burn Velocity Dispersions 


Defining shaping matrices similar to Event 2,* 


yY A Ve A 3 
M;, LT and My. SS We xe | 
37*3 3"*3 
Equations 36 and 37 become 


dt3 _ M;, O39 O21 d (Ax) 


dx = My. 83, O21 d (Ax; ) 
And again, substituting Ax, from Equation 11 yields 


03x41 
dt3 = M/,®3,0 O21 


d(Av;) 
T (eis, ) : uy 
0 
dX = My. O39 O21 ; ~ . d (Avy). 


Similar to Equations 27 and 28 for Event 2, Equations 42 and 43 can be formulated as 


dt3 = M;, 0X3 


dx; = MX,0x3. 
Returning to Equation 34, 


dxf = dx; + d(Axs). 


*M” is used here for the apogee shaping matrix instead of S” because the apogee shaping matrix for the state takes 
a downstream event, such as Event 4. 
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(39) 


(40) 
(41) 


(42) 


(43) 


(44) 
(45) 


(46) 


one form (defined by IM.) when examining the burn event itself, and a different form (defined by Sj.) when examining 


The quantity d(Ax3), however, is a function of d(Av,) and d(Av3). Therefore, 


dxt = na d( Av) a ies d(Avj,) = aren d(Avs3) (47) 
= Tox + ae | qa (Avi) + are d( Avs). (48) 
A B 


Before proceeding, it is worthwhile to understand the physical effects of d(Av;,) and d(Av3) on 
d(Ax3). B is the more obvious effect, representing the change in Ax3 due to a change in Avs. A, in 
contrast, represents the change in Ax3 due to a change in Avy, even if Av3 is zero. This is because 
a change in Av, causes a change in the inertial location of apogee at Event 3, thus changing the 
LVLH frame at Event 3. Since all burns are expressed in the LVLH frame, a change in the LVLH 
frame translates to a change in the inertial burn Av, and therefore a change to Axs. 


Returning to Equation 48, 


O(Ax3) _ O3x1 (49) 
(Avs) Tp (r3, v3 ) - Og 
Ox, O3x1 
aAr = MX,,032 021 ~ ac (50) 
ical T7(t1,V, )- tw 
and the only unknown partial now is 0(Ax3)/0x, . To obtain this partial, first note that 
O(Ax3) 0 O3x1 0 O3x1 61) 
OX; OX; T7(r3, V3 ) : U3 Av3 OX; T3 (Av3), | 
T's itself is given by 
T3 = | (hy x r3) | —hy | —fs |, (52) 


where h, = r3 X v4 is the pre-burn instantaneous angular momentum vector of the orbit. This form 
of the transformation matrix can be deduced from Figure |. Temporarily dropping the 3 subscript, 
the non-trivial derivative in Equation 51 1s 


O(T(Av),)  O(h-xr) oh- OF 
ee er ee one ay es ae 
Ox Ox oy ), Ox ( 
The (Av) components are outside the partial derivative since (Av). is not a function of x”. Re- 
calling that 


(53) 


- = ‘6 
x Sa avy (54) 
the three partial derivatives on the right-hand side of Equation 53 are given by 


Ob- ft b-b-! 
Ox- = |b7| 


I3x3 — ——_ | Ih~ x]+[rx|[v~ x] : —(rx]? (55) 
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Oe eae (Tox: ss 7) | —~[v-x] | [rx] | (56) 


OT 1 rr- , 
a (Tox: aa =) | T3x3 | 03x3 } (57) 


where b~ = h7 x r and [ax] (ais r, v, or h7) is the cross product matrix 


0 —A3 a2 
lax| = Q3 0 —Q]1 |. (58) 
—ag 1 0 
Therefore, 
0 
O(Ax3) mn x 2x0 
=| ab- Oh OF (59) 
Ox Pera He 25 Soe 25 ey 
3 axz | V32), axz | Vay), axa | V32), 
Events 4 through 7 


Events 4 through 7 work in much the same way, using the concepts and equations developed for 
Events | through 4. The result for Event 7 is 


dx7 = St bxz, (60) 


where 


O(A 
0X7 = (#0 S%, | (oa Si, )®s, ( x) Es 


O(Av,) 
O(Ax3) O(x; ) O(Axe) wae 
©-¢ 5S). |@ d(A 
(Pr _ °° Oxy Daa ee ~ (Avi) es 
(61) 
O(Axsz) O(Ax6) O(xg ) | 
@2-¢6S!. |® cy d(A 
( 7,6 i) 6,3 aha 7,6 Oxz B(Avg) (Av3) + 
O(Ax6) 
5 a) [AO 
and Sj. is defined by 
ees kf 
Dac = Per a (i *) | k = 3, 6. (62) 
k~ X;, 


The apogee shaping matrix S7 is different from MJ discussed in Event 3, thus requiring a change 
in notation. In the case of Event 3, this difference emerges from the fact that dxz (and therefore 
M ea) is not a function of x whereas dx, (and therefore Sx.) is a function of xo. In order to 
arrive at Event 4, first a time slip in t3 must be applied in order for the burn at Event 3 to occur 
at apogee, as discussed previously. If the trajectory simply stopped there, no further action would 
be required and MI a would be sufficient. To move forward from Event 3, however, dt3, must 
effectively be undone in order to maintain time coherence from Event 3 to Event 4. 
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Event k 


Equations 60 and 61 can be generalized to determine the state differential at an arbitrary Event k. 
Assume there are n burns prior to t,. The 7 burn (i = 1,...,7) occurs at t f(a) < te, where f (7) 
is the event number on which burn 2 occurs. Using the above example, with n = 3 and k = 7, the 
burns 7 = 1, 2,3 occur at events f(z) = 1, 3,6, respectively. If Event & has a constraint, then 


iO), Oy (63) 
<p, = o. OXh, (64) 
where 0x; is a lengthy expression that is omitted for brevity. The general pattern of 0x,, however, 
can be inferred from Equation 61. 
Constraint Partial Derivatives 


The analysis now turns to the LAN and SMA altitude constraint partial derivatives, which require 
the state partial derivatives developed in the previous section. The goal is to develop the partial 
derivatives of LAN and SMA altitude with respect to the burn Av magnitudes. 


LAN (A) is a function of position in the TOD frame (r,,,,) and time (¢). Temporarily dropping 
the TOD subscript for notational simplicity, the x and y components of r in the TOD frame are 


x = Pr cos @cos 6, (65) 
y= rcos¢sin9@_, (66) 


where ¢ is the geocentric latitude and 0; is the local sidereal time. #7, in turn is given by 
G7 =e, (67) 


where 6g is the Greenwich mean sidereal time and J is the longitude. Without loss of generality, 
set @ = O since A is not a function of ¢. The partial derivatives of A with respect to x and y are then 


OX y ON &£ 
eS d —=-5. 68 
Ox pe Oy rr? on 
Resuming the TOD notation, 
OX OX OX 1 
OXrop 7 Orop ON 55 r2 | ee moon : , , : ; 7 
Therefore, 
a) 1 TOD 0 
Ae = my s oy) | —Yrop “trop V9 VO OV DO —. ee (70) 
3X3 J2K 
Returning to Equation 67, 
A = 0, — 0g = O, — we (t — to) — 9G), (71) 


where Wwf 1s the magnitude of Earth’s rotation rate in rad/sec, to is an arbitrary reference time, and 
Aq, 1s the value of Og at to. Therefore, 
OX Z 


ay = WE: (72) 
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Using the chain rule, the total differential of \ at Event k is 


dA = ee: ale ai = mi, Sx, —_ WE Si, OXk, (73) 


from which the partial derivatives for OA;,/0 (Av f()) can be obtained. 


SMA altitude, in contrast, is given by a — R,. Beginning with the vis-viva equation, 


o [Ht [ 
2 69° 2a" 74) 
it is straight forward to show that 
Oa 1, le 
os mi 4 9a? | “at ad } (75) 


Since SMA altitude constraints do not need to occur at a condition such as the ascending node 
or apogee, no shaping of the state dispersions about a nominal trajectory is required. Therefore 
dx; = 0x; for SMA altitude constraints, and the total differential of a at Event k is 

Oar sh 

es = m,,,OXk; (76) 
from which the partial derivatives for Oa,/0 (Av t()) can be obtained. Note that this formula- 
tion assumes an osculating or Keplerian value for SMA. In practice however, the ISS trajectory is 
planned using a mean value for SMA, given by 


3 “| 


daz, = 


3 
a=p- rs 1 — = sin7i + sin’é cos 2u (77) 


2 
where p = a(1 — e”) is the semi-parameter and u is the argument of latitude. Equation 77 comes 
from internal TOPO documentation and is stated here without attempt to verify its pedigree. The 
partial derivative of a with respect to x 1s is given by 

Oa 7 xn Oa|Op0a Op Oe oo oe! 

—$ =-—- Mm. = — =, — — a = 

Ox “Op |0adx Oe dx Oi Ox =Ou Ox’ 
where the partial derivatives on the right-hand side of Equation 78 have been omitted for brevity. 
Similar to osculating SMA, the total differential of a at Event k is 

- Oar, T 
dap = re: = m5, OXk (79) 


from which the partial derivatives 0a, /0 (Av f(a) can be obtained. 


(78) 


The preceding developments for A, a, and a can be extended to other constraint functions as 
needed. Assume c is a generic, scalar constraint function at Event & of the form 


Ch = f (Xx, tk). (80) 
If cy, is differentiable with respect to both x and 7, then 
Oc; Oc, 
dc, = | ——S,, + —S;, | 6 81 
Ch = |e ee Bt, te | Ome (31) 


from which the partial derivatives Oc;,/0 (Av fi) can be obtained. In this expression, S,, and 
S;, are the appropriate state and time shaping matrices for cz, respectively. These matrices may be 
identity, zero, or non-trivial depending on the nature of cz. 
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REAL-WORLD ISS PROBLEMS AND RESULTS 


The above framework is now used to analyze two ISS trajectory event sequences. For each se- 
quence, the ISS trajectory is numerically integrated from event to event using FreeFlyer v7.2.1.37946 
under the following assumptions. 


e Force model: Non-spherical Earth gravity (EGM-96 7x7), Sun and Moon point mass gravity 
(DE430), and atmospheric drag (analytic density model, based on a constant temperature 
hydrostatic model). 


e Docking and undocking events are ignored, and the ISS total mass, coefficient of drag (Cd), 
and drag area are fixed to constant values. 


e All burns are assumed to have yaw = pitch = 0°. 


FreeFlyer numerically integrates the six-element state (position and velocity) along with the thirty- 
six element STM from event to event, saving off the state, STM, and other necessary results at each 
event. After reaching t+, a FreeFlyer script performs chain rule operations on the saved results in 
order to form the constraint gradients. Given the time frames spanned by ISS look-ahead trajectory 
plans—six weeks to roughly two years—obtaining the constraint gradients in this manner (via single 
propagation from fo to ty) is generally more efficient than using finite differences. 


This entire integration process can be thought of as a black box, which is called by MATLAB’s 
fmincon function to solve the constrained burn optimization problem. fmincon begins with an initial 
estimate of the Av values (discussed further below), which are the controls in the optimization 
problem. At each iteration, MATLAB sends the current Av values to FreeFlyer, which in turn 
updates the event sequence with these values, re-integrates the ISS trajectory from to to t¢, computes 
the constraint gradients at t, and returns the results to MATLAB. fmincon then uses these gradients 
along with the objective gradient in Equation 2 to update its estimate of the Av values. This process 
repeats until MATLAB converges on a locally-optimal (i.e. minimum total Av) solution, assuming 
one exists. A feasible solution—one that simply adjusts the Av values to meet the constraints—can 
also be obtained by simply setting J = 0 instead of using J from Equation 1. 


The first sequence to be analyzed consists of events taken from a TOPO look-ahead plan gen- 
erated in February, 2016. The relevant events (1.e., those having burns or constraints) are listed in 
truncated form in Table 1. In this table as well as in Table 2 for the second sequence, “sm” means 
a reboost performed using the ISS Service Module (SM) main engines, “#s” means a Soyuz event, 
and ‘“‘#p” means a Progress event. The initial Av values in Table | are taken directly from the look- 
ahead plan, while the final Av values are the optimized values. Figure 8 shows the mean SMA 
altitude plot using both initial and optimized Av values. Figures 9 and 10 show the total Av vs. 
iteration and the maximum constraint violation vs. total Av, respectively. 


In order to avoid event overlap, the second sequence consists of events from a different TOPO 
look-ahead plan generated slightly more than a year later, in March, 2017. The table and figures 
that follow mirror those used in the previous sequence. One difference in this sequence is that the 
initial Av values in Table 2 are obtained using rules of thumb that relate Av to the value of the next 
downstream constraint. These rules of thumb and their derivation have been omitted for brevity. 


Examining the first event sequence, Figure 9 is suspect since the initial estimate, taken directly 
from the actual look-ahead plan, results in an infeasible point, and the optimal total Av is in fact 
higher (8.084 m/s) than the initial total Av (7.900 m/s). This effect, however, is attributed to dif- 
ferences in force models, propagators, environment parameters (e.g. Earth’s rotation rate), event 
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Figure 8: Mean SMA Altitude Plot, February through Mid-August, 2016. The final burn and SMA 
altitude constraint “normalize” the solution so that both the initial estimate and optimal solution end 
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Figure 9: Total Av vs. Iteration 
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Figure 10: Max Constraint Violation vs. Av 
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Figure 11: Mean SMA Altitude Plot, March through Early September, 2017. Normalized solution. 
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Table 1: ISS Trajectory Events, February through Mid-August, 2016 


Initial Av Final Av Min LAN/Alt Max LAN/AIt 
Name Date Type (m/s) (m/s) (deg/km) (deg/km) 

6l1p_reboost3 01/27/2016 Burn 0.800 0.954 
6lp_reboost4 02/17/2016 Burn 0.600 0.772 
44s landing 03/02/2016 LAN 127 16.5 
63p_reboost7 08/18/2016 Burn 0.000 1.794 
SMA altitude 08/19/2016 Alt 405.5 405.5 

TOTAL 7.900 8.084 


Table 2: ISS Trajectory Events, March through Early September, 2017 


Initial Av Final Av Min LAN/Alt Max LAN/AIt 
Name Date Type (m/s) (m/s) (deg/km) (deg/km) 

sm_reboost] 03/02/2017 Burn L735 1.201 
sm_reboost2 04/05/2017 Burn 0.000 2.539 
48s_landing 04/10/2017 LAN 2 16.5 
sm_reboost8 09/05/2017 Burn 0.855 0.000 
SMA altitude 09/07/2017 Alt 405.0 405.0 

TOTAL 8.282 8.189 


triggers, etc. between this approach and the currently process for creating the trajectory plan. In 
general, direct comparison between trajectory plans generated using the approach outlined here and 
those generated using the current process is limited. This is especially so given the multi-month 
propagations that amplify these differences over time. Notwithstanding, in this case the optimizer 
takes the initial, infeasible point, and converges on an optimal solution that meets the constraints. 
The second sequence exhibits similar behavior, starting from an infeasible point and converging on 
an optimal solution subject to the constraints. However, in this case the optimal solution is at least 
not worse than the initial estimate. 


In both sequences, as noted in Figures 8 and 11, the final burn and SMA altitude constraint 
“normalize” the solution to ensure consistency. Absent this normalization, a lower total Av solution 
may be available, but the resulting SMA altitude at ty will likely also be lower. Since the ISS is 
continuously subject to atmospheric drag, this reduction in altitude, or equivalently, energy, must 
be made up at some point in the future, beyond the time horizon being analyzed. Thus any Av 
reduction achieved over [to,t 7] 1s arguably only temporary. Normalizing the solution via a final 
burn and SMA altitude constraint resolves this issue, however when coupled with the tight LAN 
constraints, it leaves little room for substantive Av savings. 


This being said, a non-normalized solution, one that does not have a final burn and SMA altitude 
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constraint, still has value for trajectory planning. In a situation where it becomes necessary to ration 
propellant, a non-normalized solution can provide non-trivial Av savings over a given period of time 
without considering longer-term altitude degradation. Such a situation 1s in fact feasible, given that 
there have been three Progress vehicle failures since 2011, each of which was carrying ISS resupply 
propellant that was obviously lost with the vehicle. While propellant rationing was not required 
in these cases, had it been required, a non-normalized minimum Av solution would provide key 
information to program managers. 


Returning to the March, 2017 event sequence, Figures 14 through 16 show this same sequence, 
but excluding the final burn and SMA altitude constraint. The optimizer in this case 1s able to 
quickly find an optimal solution that meets the constraints, and deliver almost a full 1.5 m/s of Av 
savings over the normalized solution. This represents a significant savings, given that a typical 
six-month plan requires anywhere from 6 to 10 m/s of Av. These savings, of course, will vary 
depending on several factors, such as the pedigree of the initial estimate, the number and type of 
events in the plan, and the plan duration. 


Mean SMA Altitude vs. Date 


Initial Estimate 
=== Optimal Solution 
A Launch 
7 =Burn 


V_ Landing 


Mean SMA Altitude [km] 


4' A A “A 4" A 4 4" A 4 “A 7 4" A 
RS RS RS ws RS RS KR N RS RS RS RS K RS w 
Vv v7 Vv iv v7 iV . V NY wv NY N Vv Ny 
ms AS) ms \\ \ .) no Xv No rn nS 9 No 
FS PM LK Wg é Se LP 
Date 


Figure 14: Mean SMA Altitude Plot, March through Early September, 2017. The non-normalized 
solution and its effect on the trajectory are clearly visible at t-. 


CONCLUSIONS 


This work provides a formulation and solution of the constrained burn optimization problem 
for the ISS. The formulation includes analytic derivatives for the objective function, and varia- 
tional equations using STMs for the constraint functions. These equations are used in a combined 
MATLAB-FreeFlyer framework to provide locally-minimal total Av solutions to real-world ISS 
burn planning problems. 


This framework provides several unique advantages over the current process for obtaining ISS 
burn solutions. First, it introduces automatic computation of constraint gradients with a single 
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trajectory propagation from to to t¢. These gradients can then be utilized in any gradient-based 
solver to generate feasible or optimal burn solutions for ISS burn planning. Finally, the constraint 
gradient formulation can easily be extended to other functions of the state as needed. 
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